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specific heat at constant pressure 
index of grid location in the x-direction 
index of grid location in the y-direction 
total number of grid points in the y-direction 
pressure 

radius of wall (axisymmetric case) 

Reynolds Number (per unit length) 

coefficient in front of derivatives of the advective terms m the x-directio 
coefficient in front of derivatives of the advective terms m the y-direction 
temperature 

velocity component in direction of core flow 
velocity component perpendicular to core flow 
location in the direction of the core flow 
location in the direction perpendicular to the core flow 
total viscosity (molecular and turbulent) 

core flow value at point with maximum Mach number, used to non-dimen- 
sionalize boundary layer equations 


INTRODUCTION 

At times we need to analyze the thermal behavior of systems that include both 
conduction and high speed flows. Unfortunately, most high-speed-flo^odes have 
limited conduction capabilities and most conduction codes, such as SINDA, cannot 
model high speed flows. It would be useful to interface a high-speed-flow solution and 
SINDA When interfacing a high-speed-flow solution to SINDA it may be necessary to 
include the viscous effects in the energy equations. Boundary layer effects of interest 
include heat transfer coefficients (including convection and viscous dusipafion) an 
friction coefficients. To meet this need, a fast, uncoupled, compressible, two-dimension- 
al boundary later algorithm was developed that can model flows with and without 
separation^ This algorithm was used as a subroutine with SINDA Giventhe coreflow 
properties and the wall heat flux from SINDA the boundary layer algorithm returns a 


117 


PRECEDING PAGE BLANK NCI VllMtO 


wall temperature to SINDA. SINDA and the boundary layer algorithm are iterated until 
they predict the same wall temperature. 

BOUNDARY LAYER ALGORITHM 


Boundary Layer Equations 


The forms of the boundary layer equations used in the finite difference scheme 
were the compressible, parabolized Navier-Stokes (PNS) equations by Roach, et al [1]. 
The equations were 2-dimensional, viscous, and were solved for the primitive variables. 
The y-momentum equation simply reduced to 
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The other equations were as follows: 
continuity 

d(p «) , d(p v) _ Q 

dx dy 

x momentum 


energy 
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The variables were defined as: 
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P mtl^wm 


For axisymmetric flows, (k= 1), R« was defined as: 

R o -r±ycasB 


( 10 ) 


where r was the radius of the wall, 0 was the angle of the wall, and the plus or min 
signs referred to external or internal flows respectively. For rectangular coordinates R 0 
was ignored (k=0). The coordinates x and y used in the above equations were the 
transformed x and y coordinates, that is, the * coordinate followed 
and the y coordinate was perpendicular to the wall. It should be noted that Re 
actually non-dimensional. Re had units of (length) 1 . Also, x and y (and R c ) r ® mai " e ^ 
dimensional. So, each term in each of the above boundary layer equations had units of 

(length)'. 


Solution Algorithm 

The momentum, continuity and energy equations were differenced as described by 
Kwon et al.[2] For any scalar quantity <t> : 


NM Mi-U (ID 

l&Jy x r x i - 1 

djM Vi^-1 (12) 
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± yj.ry t *-*■■ ( 13 ) 


Due to the parabolic nature of the boundary layer equations, the governing 
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equations were solved by marching from station to station in the direction of the core 
flow. Solutions at a given station were obtained by solving the boundary layer equations 
sequentially. First, the momentum equation was solved for the velocity component in 
the core flow (streamwise) direction (u). Second, continuity was solved for the velocity 
component in the y direction (v). Then, the energy equation was solved for the tempera- 
ture (T). If the station had not converged, the momentum, continuity, and energy equa- 
tions were solved again for the velocity components and temperature, using the flow 
properties from the previous iteration. Convergence at a given station was obtained 
when the streamwise velocity components at all grid locations at that station converged. 
After the calculation at the station converged, the algorithm marched to the next stream- 
wise station to solve for the boundary layer properties. This streamwise marching 
continued throughout the entire domain. 

Grid spacing perpendicular to the wall (y-direction) was based on an exponential 
function. The grid spacing was fine near the wall to better resolve the gradients at the 
wall. The grid spacing was course away from the wall where fine resolution of the 
gradients normal to the core stream were not necessary. 

The derivative of the velocity at the wall (used to calculate wall shear stress) was 
determined using a second order approximation. Values at j = l (wall), j=2, and j=3 
were used. 


~ (14) 

The above algorithm worked well for flows with weak viscous/inviscid interaction, 
since the downstream influence could be neglected. However, when the flows were 
strongly interacting, such as those with strongly adverse pressure gradients or separation, 
downstream conditions had to be considered. To account for downstream influences the 
pressure gradient was differenced as a weighted average of forward and backward differ- 
ences. This techniques was based on the method of Davis and Bamett[3]. 

Specifically the pressure gradient was differenced as 


BP 

dx 
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where c is a weighting parameter of the forward and backward differencing. This term 
was required to remove the ellipticity in the PNS equations in strongly interacting flows. 
The quantity £ determined what fraction of the forward difference of the pressure 
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gradient can be included so that the equations remain non-elhpud If the 
supersonic at a given j location, then only backward differencing was used for M*e p 
sure gradient (e = 1). If the flow was subsonic then at a given j location then th 
following expression for c was used, 


e= 


yM 2 


l+iy-W 1 


(16) 


Boundary Layer Separation 

Haring of the advective term parallel to the core flow was used when the flow 
was separated. This was handled by taking the absolute value rffte c«fBaent in front 
of the derivative of the advective term (RUB) and multiplying it by .1 to make it smaller. 
Again this operation was only performed when the flow was separated. 

Boundary Layer Turbulence Model 

A modified Baldwin-Lomax model was used to account for turbulence. This 
model was a zero equation, eddy viscosity model. This model w^ d 

turbulence models, such as the k-e two equation model. A modified model of Vubal 
Knight was used to handle the effects of separation. However the model of Visbal and 
Knight required further modification, the most important of which was the modification 
of the Baldwin-Lomax parameter C^. C,, needed to be a function of both i core : flow 
Mach number and core flow pressure gradient. This modification was extremely 
important in matching numerical results and experimental data. A detailed description 
of the turbulence model used is given in Sakowski, et. al [4J. 

Boundary Layer/Core Flow Interface 

The boundary layer algorithm serves as a link between the conduction program 
(SINDA) and an inviscid core flow program. The interface with SINDA will be is- 
cussed in a later section. In this section we will look at the interface of the boundary 
layer algorithm with a core flow program. 

When the boundary layer algorithm was interfaced with a core flow algorithm, it 
was necessary that the boundary layer properties smoothly approached the core flow 
values. That is, when the derivatives perpendicular to the waUwere zero, the core flow 
values had to be a solution to the boundary layer equations. The complicated part was 
ma tching the boundary layer and core flow, since the boundary layer and core flow algo- 
rithms were probably not differenced in the same way. If the differential equabons were 
solved exactly there would not be problem, but they were not solved exactly. What the 
core flow algorithm predicted as a solution, was not exactly what the boundaiy layer 
algorithm predicted as a solution as the y-derivatives went to zero (far from the wall). 
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The difference was usually fairly small (2% or so), but this small -difference could have a 
big effect on the integral performed to calculate the displacement thickness. An 
adjustment in the way the pressure gradient term was calculated in the boundary layer 
algorithm forced the boundary layer properties to smoothly approach the core values. 
Without this adjustment, the displacement thickness, predicted by the boundary layer 
program had large errors. The adjustment of the pressure gradient was performed by 
solving the finite differenced momentum equations for dP/dx when all y-derivatives were 
zero. 


dP 

ax 



(17) 


The corrected value of the pressure derivatives was calculated with the edge 
values from the core flow algorithm using the same differencing scheme used in the 
boundary layer algorithm. In this way the boundary layer algorithm approached the core 
flow values as the y-derivatives approached zero. 

Another consideration for interfacing a core flow algorithm and the boundary 
layer code was stability. At times the boundary layer algorithm has a stability problem. 
This problem tended to initiate near the edge of the boundary layer. From one iteration 
to the next the values near the core flow sometimes fluctuated between less than and 
greater than the core flow value. Sometimes these fluctuations died out and the program 
converged. However, other times the oscillations grew, causing the calculations to 
diverge. To solve this problem flaring was used. RVB was part of the advective terms 
in the y-direction. RVB was the coefficient in front of the du/dy term in the finite 
differenced x-momentum equation, and the dT/dy term in the finite differenced energy 
equation. These were the advective terms in the y-direction. Without flaring RVB was 
simply pv. With flaring RVB was changed as follows: 

RVB=Kj p v | < 18 > 

K,’s for momentum equation 

if u<u t if u>u e 


if — >0 

ay 




if ^<o 
dy 





u 
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K,’s for energy equation 


if T<T t if T>T t 



The motivation for the above flaring was to make the core flow value a numeri- 
cally stable solution in the boundary layer algorithm far from the wall. This flanng was 
found to be very important to help the stability of the algorithm, particularly when there 
was an adverse pressure gradient, separation, bleeds, or bypasses. For more detail on 
interfacing the boundary layer algorithm of Roach, et al. with a core flow, refer to 
Darling, et al. [5]. 


CORE FLOW INPUT TO BOUNDARY LAYER ALGORITHM 

It was mentioned previously that the boundary layer algorithm required as an 
input, a core flow. The inviscid core flow variables needed by the boundary layer 
algorithm were the Mach number, temperature, and pressure. For flows where the 
interaction between the core flow and boundary layer were negligible, the interface was 
very simple and easily implemented. All that was required was an mput file by the name 
LBL.DAT. This file contains seven namelists described as follows: 

NAMELIST/GEOM/X: Grid point locations that run parallel to the centerline 

of the wall 


NAMELIST/AREAX/AREA: 

NAMELIST/TEMP/T1: 

NAMELIST/PRES/PI: 


Surface area of the wall that corresponds to 
a respective X grid point location 

Core flow temperature that corresponds to a 
respective X grid point location 

Core flow that corresponds to a 
respective X grid point location 
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NAMELIST/MACH/AM: Core flow Mach number that corresponds to a 

respective X grid point location 

NAMELIST/RADIUS/RADY: Radius of the wall measured from the wall centerline 

and perpendicular to its respective X grid point 
location 


NAMELIST/BL/ 

jjj: Number of grid points perpendicular to the wall 

DUMX: Convergence criteria on the u velocity component of 

the x momentum equation 

II: Number of grid points in the streamwise direction 

TURB: Flag to signal use of the Baldwin-Lomax turbulence 

model: 

TURB = .' TRUE. - turbulence on 
TURB = .FALSE. - turbulence off 

ENGU: Flag to signal use of units: 


ENGU 

T1 

PI 

X 

AREA 

ENG UNITS - .TRUE. 

•R 

lbf/ft 2 

ft 

ft 2 

SI UNITS - .FALSE. 

•K 

N/m 2 

m 

m 2 


AXI: Flag to signal axi symmetric flow: 

AXI = .TRUE. - axisymmetric flow 
AXI = . FALSE. - 2-D flow 

EXT: Flag to signal external flow: 

EXT= .TRUE. - external flow 
EXT =. FALSE. - internal flow 

MYES: Flag to print to: 

Mach.out: Prints non-dimensional velocity, 

temperature, density, and pressure 
profiles for each x grid location and for 
the first 10 SINDA iterations - after 10 
SINDA iterations, the profiles are 
printed eveiy 20 SINDA iterations 
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Each x grid location also gves the 
following parameters: 

IT: number of boundary layer 

iterations 

M: core flow mach number 

Cf: skin friction coefficient 

Ybl: distance of the last y grid point 

from the wall 


C2.0UT: Prints the following boundary 

layer parameters: 

IT: number of boundaiy layer 

iterations 

Th: momentum thickness 

Cf: skin friction coefficient 

Disp: displacement thickness 

for the first 10 SINDA iterations - after 

10 SINDA iterations, the profiles are 

printed every 20 SINDA iterations 


A sample LBL.dat file is provided in the APPENDIX. 

For flows where there was a strong interaction between the core flow and 
boundary layer the simple input data file of constant inviscid core flow values would 
change complexion. The file would become a core flow algorithm which continuously 
updates the rare flow variables to account for the boundary layer interaction, sue 
flows with shock waves. Such algorithms are not discussed in this paper. 


SINDA/BOUNDARY LAYER INTERFACE 

The nodes in SINDA that represent the wall surface nodes MUST be 
boundary nodes. This is done to obtain heat rates on the nodes that can be sent to the 
boundary layer algorithm. For a steady state solution the heat rate on the surface 5 node 
will be^-ero if the nodes are defined as arithmetic or diffusion nodes. This “ 

SINDA effectively "sees" an insulated surface. In actuality the surface is ’ 

because of the presence of the boundaiy layer flowing over it. For each SINDA 
iteration, the heat rates from the SINDA nodes are passed to the bounda^^yer where 
new wall temperatures are determined. The boundaiy layer algorithm willcontinue 
iterating until it has a converged solution. After the boundaiy . , 

“ed the bcuodary layer wall temperatures become the new SINDA Mr-g. 
temperatures. Thus the SINDA boundary node temperatures are updated eveiy SINDA 
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iteration. The program iterates between SINDA and the boundary layer algorithm until 
the boundary layer wall temperatures match the SINDA boundary node temperatures. 


SINDA INPUT TO BOUNDARY LAYER ALGORITHM 


In HEADER CARRAY DATA, reserve the following variables: 

1-LBL.DAT 
2= MACH. OUT 
3=C2.0UT 

Also, the common SINDA/FORTRAN variable BTEST should not be used. In HEAD- 
ER VARIABLES 1, make the call to the subroutine INTERFACE as follows: 

CALL INTERFACE(A,B,C,D,E,ABSZRO,BTEST,UCAl,UCA2,UCA3) 

where A = Submodel name in quotes where the surface boundary nodes are located, 
and no more than 8 characters long. 

B = SINDA node number of the surface boundary node located at the last x 
grid location of the boundary layer (integer). 

C = SINDA node number of the surface boundary node located at the first x 
grid location of the boundary layer (integer). 

D = Number by which the SINDA boundary nodes are incremented (integer). 

E = Units used by SINDA. ’ENG* for English units. ’SI’ for SI units. 

The remaining arguments should be left as they are. ABSZRO is the SINDA variable 
for absolute temperature defined in the HEADER OPTIONS DATA BLOCK. BTEST 
is a counter for the number of SINDA iterations. UCA1, UCA2, and UCA3 are the 
SINDA variables for the input and output file names defined in the HEADER CAR- 
RAY DATA BLOCK. 


CONCLUSIONS 

A fast steady, compressible, turbulent boundary layer algorithm that can be used 
to model separated flows has been written as a subroutine for SINDA Results from the 
boundary layer algorithm compared well with experimental pressure distributions when 
the boundary layer was interactive with the core flow, Darling, et al. [5] and Roach, et al 
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[11. In addition, the friction coefficients and momentum thicknesses predicted by the 
boundary layer code compared well with experimental data. Roach, et al. [1]. 
EkddwhvLomax turbulence model was used following the modifications of 81 

[41 also matched experimental data fairly well. Currently comparisons are bcl ^g ™ ad 
with experimental data validate the heat transfer predictions of the boundary layer 
algorithm. The boundary layer algorithm was found to cower* ' 
simple SINDA model with 25 x-grid locations was tested, see the APPENDIX for the 
SINDA input. The model converged in 95 SINDA iterations. 
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APPENDIX 


SGEOM X* 


225528E-2, .247675E-1, .515367E-1, .778460E-1, .107633, .128037, .151619, 
174448, .203611, .228309, .231856, .235408, .238954, .242498, 


.246048, .249659, 

.25352, .258018, .262899, .267638, .272238,-276733, 

.281173, .285603, .290036, .294476, .298924, .303379, .307842, .312309, 
.316781, .321258, .325737, .330218, 

.334703,. 343673,. 352644,. 366070 SEND 

fcAREAX AREA«1.0 # 1.0,1.0,l*0,1.0 # 1.0 # 1.0 # 1.0 f 1.0, 1.0, 1.0 ,1.0, 1.0,1 •0,1. 0,1.0, 

1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 

1.0, 1.0, 1.0, 1.0, 1.0, 1.0 SEND 

STEMP Tl*1982 . 9 , 1983 .3,1983.5, 1983 . 6 , 1983 . 5 , 1983 . 3 , 1982 . 7 , 1981 . 9 , 1980 . 3 , 

1975.5.1973.5.1970.8.1966. . 1960.4.1945.2.1917.8.1857.5.1776.7, 
1707.9,1605.9,1526.4,1468.8,1385.8, 

1325. 1.1283.. 1253. 3. 1228. 5. 1207. 9, 

1188.9.1171. . 1155.3.1142.7.1131.3.1122.4, 


1116.5,1106.3,1101.8,1098.6 SEND 

fiPRES PI*. 172115E7 , . 172253E7 , . 17236E7 , . 172399E7, . 172346E7 
. 17199E7 , . 17161 5E7, . 170832E7 , . 168388E7, . 167384B7 


172235E7 , 

165993E7, .163667E7, 


. 160348E7 , . 1S3294E7 , .141769E7, .122549E7, 

.101028E7, .810816E6, .613865E6, .492139E6, .407487B6, 

.315702E6, . 257376E6, .222281E6, .20018E6, . 183118E6, .169834E6, 

. 1S8537E6, . 148288E6, . 139641E6, . 133107E6, . 127322B6, . 122975E6, 


. 120144E6, . 115385E6, 

•113278E6, .111714E6 SEND 

&MACH AM*. 389E-1, .393E-1, .411E-1, .442E-1, .492E-1, .541E-1, .619E-1, .73E-1, 

.105, .192, .217, .25, .29, .354, .442, .561, .725, .919,1.13,1.3,1.48,1.6, 
1.77,1.89,1.98,2.01,2.1,2.14,2.19,2.22,2.25,2.29,2.31,2.32,2.34, 

2*36 2*37 2*374 SEND 

SRADIUSY RADY*1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1. 0,1.0, 

1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 

1.0, 1.0, 1.0, 1.0, 1.0, 1.0 SEND 
SBL JJJ-100, DUMX1*. 001, II«25,TURB*. TRUE. , 

ENGU-. FALSE. ,AXI*.TRUE. , EXT*. TRUE. , MYES*.TRUE. SEND 
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title R boundary layer code interface with sinda 

MODEL - BL1 
OUTPUT « BL2.0UT 
USER1 * BL1.USR 
USER2 « BL2.USR 
HEADER CONTROL DATA, GLOBAL 
NLOOPS * 4000 
ABSZRO * -460.0 
UID ■ ENG 
ARLXCA - .01 
DRLXCA « .01 
EBALSA « .01 

HEADER USER DATA, GLOBAL 
C 

HEADER CARRAY DATA, BL 


JDEFAULT VALUE 
JDEFAULT VALUE 
SDEFAULT VALUE 


l-NASASPFSD: 

2«T0ISK$DIR: 

3*TDISK$DIR: 


[AMBER. SIN0A]LBL. DAT 
‘AMBER] MACH. OUT 
‘ AMBER] C2. OUT 


HEADER USER DATA, BL 

C»»»RESERVE BTEST FOR USE IN BOUNDARY LAYER 
HEADER ARRAY DATA, BL 
C 


81. 

,0.1139 

261. 

,0.1230 

621. 

,0.1330 

981. 

,0.1390 

1341. 

,0.1459 

81. 

, 8.61 

261. 

, 9.59 

621. 

,11.44 

981. 

,13.06 

1341. 

,14.68 


$ SPECIFIC HEAT VS. TEMPERATURE 
$ UNITS: BTU/LBM/DEG. F 

$ AISI 304 S.S. 


$ THERMAL CONDUCTIVITY VS. TEMPERATURE 
S UNITS: BTU/HR/FT/DEG. F 

S AISI 304 S.S. 


C 

HEADER NODE DATA, BL 


r******AISI 304 S.S. /DENSITY-493 LBM/FT**3 
£****«*G£N£RATE 25 DIFFUSION NODES TO REPRESENT THE WALL! 

SIM 8801,25,1,70. ,A1,5. 

c ******GENERATE SURFACE NODES! !!!!!! !»>CONNECT TO THE BOUNDARY LAYER CODE 
GEN -1,25,1,3000. ,0.0 

c ****** GENE RATE OTHER BOUNDARY NODES TO SIMULATE THE EFFECT OF A SIMPLE 
e ******ACTI VE COOLING SYSTEM WHOSE EFFECTIVE TEMPERATURE IS 100 DEG FI!! 


GEN -1101,25,1,100.0,0.0 


HEADER CONDUCTOR DATA, BL 


r********rDPATf CONDUCTORS IN THE WALL ALONG THE "X AXIS" _ ci nu 

cist THE •" AXIS- FOLLOWS THE DIRECTION OF THE CORE FLOW 

SIM 801, 24, 1,8801,1, 8802, 1,A2, 6. 
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r****«*** CREATE CONDUCTORS THAT CONNECT THE WALL DIFFUSION NODES , 
c ****«***TO THE BOUNDARY NODES WHICH INTERFACE WITH THE BOUNDARY LAYER 
SIM 8801,25,1, 8801,1,1, 1.A2..1 

C ******** C REATE CONDUCTORS THAT CONNECT THE WALL DIFFUSION NODES 
C ********T0 THE BOUNDARY NODES THAT SIMULATE ACTIVE COOLING 
SIM 88801, 25, 1,8801,1, 1101, l,A2,.l 

HEADER VARIABLES 1, BL 

C ***** THIS ls the SUBROUTINE THE DOES IT ALL!!!!!!! 

CALL INTERFACE ('BL ' ,25,1,1, 'ENG' ,ABSZRO, 

& BTEST , UCA1 , UCA2 , UCA3 , ) 

HEADER OPERATIONS DATA 
C 

BUILD BL1.BL 
C 

CALL STDSTL 


# 

« 


# 
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